#include "SquareMatrix.hpp"

SquareMatrix::SquareMatrix( unsigned int n ) :
	Matrix( n, n )
{
}

SquareMatrix::SquareMatrix( const Matrix& m ) :
	Matrix( m )
{
	if( !m.isSquare() )
		throw std::logic_error( (char*)"La matriz no es cuadrada" );
}

SquareMatrix::~SquareMatrix()
{
}

// traspone la matriz
void SquareMatrix::transpose()
{
	__MATRIX_EACHCOL(j)
		for( unsigned int i = 0 ; i <= j ; i++ )
		{
			double a_ij = _matrix[i][j];
			_matrix[i][j] = _matrix[j][i];
			_matrix[j][i] = a_ij;
		}
}

// calcula el determinante de la matriz
double SquareMatrix::det() const
{
	double result = 0;

	// caso 1x1
	if( _cols == 1 )
		result = _matrix[0][0];
	// caso 2x2
	else if( _cols == 2 )
		result = _matrix[0][0]*_matrix[1][1] - _matrix[0][1]*_matrix[1][0];
	// caso 3x3
	else if( _cols == 3 )
		result = 
			_matrix[0][0]*_matrix[1][1]*_matrix[2][2] +
			_matrix[0][1]*_matrix[1][2]*_matrix[2][0] +
			_matrix[0][2]*_matrix[1][0]*_matrix[2][1] -
			_matrix[0][0]*_matrix[1][2]*_matrix[2][1] -
			_matrix[0][1]*_matrix[1][0]*_matrix[2][2] -
			_matrix[0][2]*_matrix[1][1]*_matrix[2][0];
	// caso n x n
	else
	{
		unsigned int i = 0;
		bool add = true;
		__MATRIX_EACHCOL(j)
		{
			double val = _matrix[i][j] * minor(i,j).det();
			if( add )
				result += val;
			else
				result -= val;
			add = !add;
		}
	}

	return result;
}

SquareMatrix SquareMatrix::transposed() const
{
	return Matrix::transposed();
}

SquareMatrix SquareMatrix::minor(unsigned int i0, unsigned int j0) const
{
	return Matrix::minor(i0,j0);
}

SquareMatrix SquareMatrix::LUfactorization()
{
	SquareMatrix L( _rows );

	for( unsigned int j = 0 ; j < _cols-1 ; j++ )
	{
		// pivoteo parcial
		unsigned int k = j;

		// busca la fila con mayor coef. (en módulo) en la columna j
		for( unsigned int i = j+1 ; i < _rows ; i++ )
			if( std::abs(_matrix[i][j]) > std::abs(_matrix[k][j]) ) 
				k = i;

		// si había una, intercambia con la actual
		if( k > j )
		{
			swapRows( j, k );
			L.swapRows( j, k );
		}

		// triangula la columna si no está triangulada
		if( std::abs(_matrix[j][j]) > 1.0e-14 )
			for( unsigned int i = j+1 ; i < _rows ; i++ )
			{
				double k = _matrix[i][j] / _matrix[j][j];
				L._matrix[i][j] = k;
				substractRows( i, j, k );
				_matrix[i][j] = 0;
			}
		else
			_matrix[j][j] = 0;
	}

	return L+SquareMatrix::identity(_rows);
}

void SquareMatrix::triangulate()
{
	for( unsigned int j = 0 ; j < _cols-1 ; j++ )
	{
		// pivoteo parcial
		unsigned int k = j;

		// busca la fila con mayor coef. (en módulo) en la columna j
		for( unsigned int i = j+1 ; i < _rows ; i++ )
			if( std::abs(_matrix[i][j]) > std::abs(_matrix[k][j]) ) 
				k = i;

		// si había una, intercambia con la actual
		if( k > j )
			swapRows( j, k );

		// triangula la columna si no está triangulada
		if( std::abs(_matrix[j][j]) > 1.0e-14 )
			for( unsigned int i = j+1 ; i < _rows ; i++ )
			{
				double k = _matrix[i][j] / _matrix[j][j];
				substractRows( i, j, k );
				_matrix[i][j] = 0;
			}
		else
			_matrix[j][j] = 0;
	}
}

SquareMatrix SquareMatrix::identity( unsigned int n )
{
	SquareMatrix id( n );
	for( unsigned int i = 0 ; i < n ; i++ )
		id._matrix[i][i] = 1;
	return id;
}

SquareMatrix SquareMatrix::Hilbert( unsigned int n,unsigned int cte )
{
	SquareMatrix H(n);
	__MATRIX_EACH_M(H,i,j) H._matrix[i][j] = 1.0f/(((i+1)+(j+1)-1)+cte);
	return H;
}

SquareMatrix SquareMatrix::ShootMatrix(unsigned int n, const Vector& y,const Vector& x,unsigned int k)
{
	SquareMatrix A = SquareMatrix::Hilbert( n, k );

	// elijo la columna j0, tal que x_j0 es distinto de cero
	unsigned int j0 = 1;
	while( x[j0] == 0 && j0 <= n ) j0++;
		using namespace std;

	double tmp = 0;
	for( unsigned int i = 1 ; i <= n ; i++ )
	{
		tmp = y[i];
		for( unsigned int j = 1 ; j <= n ; j++ )
		{
			if( j != j0 )
				tmp -= A(i,j) * x[j];
		}
		tmp /= x[j0];
		A(i,j0) = tmp;
	}

	return A;
}
